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ABSTRACT 

A systematic search for multicomponent crystal structures is carried out for five different ternary 
systems of nuclei in a polarizable background of electrons, representative of accreted neutron star 
crusts and some white dwarfs. Candidate structures are “bred” by a genetic algorithm, and optimized 
at constant pressure under the assumption of linear response (Thomas-Fermi) charge screening. Sub¬ 
sequent phase equilibria calculations reveal eight distinct crystal structures in the T = 0 bulk phase 
diagrams, five of which are complicated mult inary structures not before predicted in the context of 
compact object astrophysics. Frequent instances of geometrically similar but compositionally distinct 
phases give insight into structural preferences of systems with pairwise Yukawa interactions, including 
and extending to the regime of low density colloidal suspensions made in a laboratory. As an applica¬ 
tion of these main results, we self-consistently couple the phase stability problem to the equations for a 
self-gravitating, hydrostatically stable white dwarf, with fixed overall composition. To our knowledge, 
this is the first attempt to incorporate complex multinary phases into the equilibrium phase layering 
diagram and mass-radius-composition dependence, both of which are reported for He-C-0 and C-O-Ne 
white dwarfs. Finite thickness interfacial phases (“interphases”) show up at the boundaries between 
single-component bee crystalline regions, some of which have lower lattice symmetry than cubic. A 
second application - quasi-static settling of heavy nuclei in white dwarfs - builds on our equilibrium 
phase layering method. Tests of this nonequilibrium method reveal extra phases which play the role 
of transient host phases for the settling species. 

Subject headings: dense matter - methods: numerical - plasmas - stars: neutron - white dwarfs 


1. INTRODUCTION 

When an impure white dwarf (WD) or neutron star 
crust (NSC) is slowly cooled from above its melting tem¬ 
perature, one expects the extra compositional degrees of 
freedom are taken advantage of to form crystals which are 
more efficiently packed than phase-separated bee lattices. 
Indeed, several investigators have considered non-Bravais 
and multicomponent lattices as the possible ground state 
of a strophysical c ompact objects. One of the earliest 
was Dyson (1971), who suggested a rock salt structure 
of be and He nucle i migh t be stable. More recently, 
Kozhberov & Baiko (2012) studied cesium-chloride and 


magnesiu m-di boride structures within the Coulomb crys¬ 
tal model. |Kobyakov & Pethick (2014) have argued that 
the ground state structure above neutron drip density 
may be similar to that of the displacive ferroelectric 
BaTiOs, due to the symmetry-lowering effect of inter¬ 
stitial neutrons on a bee lattice of nuclei. A related line 
of inquiry concerns the freezing of multicomponent ion 
plasni as from the liquid state. See Medin fc Cumming| 
(2010) for a semi-analytic calculation and references to 
earlier numerical methods. One such method - clas¬ 
sical molecular dynamics - has been used exte nsively 


et al. d2007|) comi 

position (Horowitz et al. 

2007 

2009: 

Horowitz Berry 

2009 

). T'he latter of these works tea- r 


tures a 14-component, 
was annealed for ^10^ phonon cycles below the melt- 
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ing temperature. A dominantly Se (Z=34) bee lattice 
was formed, with small-Z nuclei occupying interstitial 
positions and larger-Z nuclei acting as substitutional im¬ 
purities. In addition, there was a tendency for small-Z 
nuclei to cluster together, forming an effective large-Z 
particle. In a different simulation where a nnealing was 
again ca rried out for ^10^ phonon cycles (Horowitz et 


al. 2009), phase-separated regions (microcrystals) formed 


in the solid phase. One phase was depleted in small-Z 
nuclei, while another was enriched. 

Simulated annealing is an excellent means for directly 
modeling the dynamics of crystalline systems, but it of¬ 
ten cannot access the very long timescales associated 
with the nucleation and growth of complex multicom¬ 
ponent crystal phases, due to the exponentially slow dy¬ 
namics of surmounting reaction barriers against the com¬ 
plex cooperative rearrangements needed to form such 
crystals. For example, terrestrial carbon steels, which 
typically have only 2-3 alloying elements, must be an¬ 
nealed for a minimum of ^10^^ p honon cycles (^ 10 
seconds) to find their ground state (|Atlas of Isothermal 


Transformation and Cooling Transformation Diagrams 


1977). Alt ernative methods includ ing random structure 
searching dPickard fc Needs] 2011), particle swarm opti¬ 
mization dWang et al.||2UlU), and genetic / evolutionary 
search techniques dOganov ^ Glass] |2QQ6| |Abraham &| 
Probert||2008| |Wu et al.||2U14|) have been applied with 
great success to this “crystal structure problem,” but 
have not yet been applied at the extreme conditions of 
compact astrophysical objects. When coupled with an 
appropriate description of the (fully pressure-ionized) mi¬ 
crophysics, such methods could provide a means to efh- 
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ciently search for new crystal structures in multicompo¬ 
nent WDs and NSCs, complementing the existing simu¬ 
lated annealing work. 

The existence of lower-symmetry (i.e. non-cubic) 
and/or multinary phases within WDs and NSCs could 
have several astrophysical implications. Most astrophys- 
ical calculations assume the material is a bee polycrystal 
with grain sizes small compared to the other macroscopic 
physical scales in the problem. Therefore, for example, 
the rank-four elastic tensor is averaged and smoothed to 
produce a scalar shear modulus relating the strain re¬ 
sponse to an applied s tress (one popular averagi ng pro¬ 
cedure is described by Ogata & Ichimaru ( 1990|)). The 
possibility of multiple, complicated lattice structures, 
and preferential alignment with e.g. the local magnetic 
field, would necessitate computing the full elastic tensor. 
Anisotropies, soft phonon modes, and el astic instabilities 
such a s the incipient ones described in [Engstrom et al. 
(2015) could have significant effects on elasticity- related 


fc Pons |2Qll), related quasi-periodic oscillations ( Israel 


astrophysical observables such as magnetar flares ( Perna 


& Haensel 


etalj2Q05), and possibly some pulsar glitches (C' 


lamel 


2008). It could also significantly affect the 


future observability of gravitational-wave emission, both 
i n the context of magnetar flares an d continuous waves 
(Johnson-McDaniel & Owen 2013). Grain/phase do- 
main boundaries would lead to preferred stress-failure 
locations, and on large scales might affect dissipation of 
modes involv i ng the crust such as torsion or shear modes 
(Israel et al. (2005), used to explain quasi-periodic os- 
cillations after magnetar flares) and r-modes (similar to 
the “crust freezing” scenario in|Lindblom et al. (2000)). 
These again would have implications both for electro¬ 
magnetic and gravitational wave observations. Another 
kind of implication has to do with the composition of WD 
debris disks and planetary systems, inferred fro m metal 


abundances in the accret ing WD’s atmosphere (Barber 
et al.|2012[ |Raflkov|2011|. Entering into this calculation 
is the settling rate of the nigh-Z met a! 


settling rate of the high-Z metals. In principle, this 
rate depends on the buoyancy of the settling species’ host 
phase(s) as well as the microphysics involved in ordinary 
grain growth processes, na mely interfacial ene rgies and 
grain boundary mobilities (Krill & Chen||2002). 

In this work we carry out a systematic search for the 
ground state crystal structure of three-component sys¬ 
tems at conditions relevant to WDs and NSCs. The 
main goals are I) to identify possible new phases through 
global search of the multicomponent crystal structure 
phase diagram, and 2) to determine in what contexts 
those phases might appear in WDs and NSCs through 
sample applications to layering stability. To the first 
end, we employ a popular genetic search algorithm. The 
lowest enthalpy structures found by the genetic search 
are included in bulk phase diagram calculations, which 
reveal five new complicated binary and ternary crystal 
structures, four having sub-cubic lattice symmetry. To 
the second end, we demonstrate a self-consistent coupling 
of the phase stability calculation with the basic equations 
of a self-gravitating, hydrostatically stable white dwarf. 
Several compositional instances of the newly found bi¬ 
nary phases show up as finite thickness “interphases” 
between pure bee strata in cold, He-C-0 and C-O-Ne 
white dwarfs. Additional binary phases make a tran¬ 
sient appearance in nonequilibrium settling calculations. 


as host phases for the settling species. 

2. BULK PHASE DIAGRAM CALCULATION 

This section describes a global search of composi¬ 
tion and structure, using five ternary systems of nuclei 
thought to be relatively prevalent in WD or NSC mat¬ 
ter, and covering a range of distinct “crystal chemistries.” 
The starting point is an effective Hamiltonian for com¬ 
pletely pressure-ionized matter. We work within linear 
respon se th eory - see Se ction V of |Pollock fc Hansen 
(1973), and Baiko ([2002), for example! In this frame¬ 
work, a system of point nuclei (with charges and 
static positions r^) immersed in a polarizable, charge 
compensating background of electrons has kinetic plus 
electrostatic potential energy 


V 




(27r)3 
d^k 47r 




/c^e(k) 

^ -1 


(1) 


e(k) 

d^k 47re**'''‘ 
(27r)^ k‘^e(k) 


(compare 

immediate. 


Baiko (2002) Equations 1-3). It is not 


y obvious that the above Hamiltonian in¬ 


cludes the leading order correction to the kinetic en¬ 
ergy To of the uniform electron gas (it does). This 
can be seen by expanding the kinetic energy in pow¬ 
ers of the density nonuniformity correction: T = 

To + ^ / d?r (fir' 6ne{r)G{T — r') Jne(r') + ... and keep¬ 
ing only the first two terms such that a total energy 
minimization identifies — G(k)“^ as the static response 
function of the uniform gas. Equation also contains 
all Coulomb interactions except for the infinite nuclear 
self energies. With a choice of the simple Thomas-Eermi 
dielectric function eTF(k) = 1 the integrals are 

standard ones and the Hamiltonian reduces to 


Etf — To-\-—< 


Att 
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^ 7^3 


E: - r,- 


( 2 ) 


In performing structural optimizations, one converges 
the energy by working with a supercell of volume W and 
N — 1 periodic copies thereof. If R is a primitive lattice 
vector (supercell translation vector) and p, q index the 
basis, the total energy per supercell is written 


N 




P, <7 




( 3 ) 


where ^pq = |R + rp — r^j and the prime on the last sum 
indicates that terms with ^pq = 0 are excluded. Eor this 
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work we use kinetic energy density 
1 nieC^ rx^(l + 2x‘^) 


tq-- 


Stt^ 




and screening length 


In 


(-1) 


■ meC Ue, 


( 4 ) 


1 ^ KM 

2x\ a ' 


( 5 ) 


corresponding to the relativistic, degenerate gas. Here 
Ae = hjm^c is the reduced Compton wavelength, x = 
PFlrrieC = Ae(37r^ne)^/^, (3 = x/^Jl + and a « 1/137 
is the fine structure constant. The Thomas-Fermi de¬ 
scription breaks down when the screening length local¬ 
izes electrons to within their Compton wavelength; this 
occurs for x ^ 10 {p ^ 10^ g/cc). Approaching this 
extreme relativistic limit, the ratio of screening length 
to Wigner-Seitz radius Vg = {3{Z) tends to a 
constant (for a bcc lattice of Fe nuclei, the constant is 
1.82), so we might anticipate that phase boundaries be¬ 
come stationary in this scale-invariant limit. Saturation 
of kQ^/vg at this small value is indicative of the over¬ 
screening predicted by the Thomas-Fermi model. Within 
these constraints however, the model has the advantage 
of being both reasonably accurate and computationally 
efficient, readily incorporated into a global search of crys¬ 
tal structure and composition. 

Structural optimizations using Equations [3}|^ are con¬ 
veniently carried out at constant volume - one need only 
minimize a pairwise sum of effective Yukawa interactions, 
which converges rapidly in real space for k^^/vg ^ 1. 
The simplicity of this method is conducive to a large 
basis, useful for studying interfac es, defects and fail¬ 


ure m echanisms (see for example, Horowitz & Kadau 


( 2QQ9| )). Because it will prove useful for the phase layer- 


ing calculation described later, we choose instead to per¬ 
form structural optimizations at constant external pres¬ 
sure P and minimize the enthalpy. This can be accom¬ 
plished using a modified ve rsion of the Genera l Utility 
Lattice Program (GULP) (Gale & Rohl 2003), where 
the Yukawa interaction (available as a special case of the 
“general” pair potential) is given the capability to handle 
a Uc-dependent screening length. We modify GULP’S en¬ 
thalpy per unit cell to hr f = Etf/Y+PW, and the first 
strain derivatives of the enthalpy (sufficient for steepest 
descents and conjugate gradient methods) are accord¬ 
ingly modified to 


SHtf 

dea/3 


(P - P)V M ~ z Z 

[P Po)Fc+ q7.2t/ / , ZpZq 


3klVc 


(6) 
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e‘^ko{l -h (3‘^) 
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p, q 


^ — ko^p 


R,p, q 


M/3 


e2 ^ 


-tE 


R,p, q 


^pq 




where Pq = jjdn^ is the kinetic pressure of the 

uniform electron gas. Derivatives of Htf with respect to 
gulp’s remaining degrees of freedom (fractional basis 
coordinates) are not affected by ko M ko{Vc). 


We carry out ground-state structure searches using 
the evolution ary crystal structur e prediction software 
XtalOpt rS.Q ^Lonie fc Znrek]2Qll ), together with GULP 
optimization^ It is assumed the ground state is a poly¬ 
crystalline mixture of stoichiometric compounds. We 
don’t consider solution (alloy) phases, and we consider 
only a subset of possible stoichiometries. For a given 
ternary system of nuclei A-B-C, a constant pressure 
search is performed at P = 10^^ and 10^^ GPa for 
each of the nominally 125 stoichiometries where 

n^m^i = 0...4. Search cells with small number of 
particles n P m P i are removed from the search pro¬ 
gram if they are sub multiples of larger cells, thus there 
are 98 searches per ternary system, per pressure. Each 
of these searches is run out to at least 480 optimized, 
genetically-operated-on structures, except in the case of 
single-component searches, which are run out to at least 
80 optimized structures (the first 20 seed structures are 
randomly generated). Lattice sums are done in real 
space with cutoff 20/c^^ and are expected to be con¬ 
verged to 7-8 digits. This level of accuracy is impor¬ 
tant, as we find enthalpies of competing structures can 
be the same out to 6 digits. Default XtalOpt search 
parameters are used throughout, and following the sug- 
gestions put forth in the XtalOpt implementation paper 
( Louie fc Zurek||2011 ), we benchmark the search param- 
eters by constructing Hartke plots for several relativistic 
screened-Coulomb systems (see Eigurej^. Hartke plots 
gauge the performance of a genetic search and help es¬ 
tablish a stopping criterion. Our choice of search dura¬ 
tion, previously mentioned, is in part motivated by the 
“Hartke lifetimes” found in these tests. 

The lowest enthalpy structure found in each search 
is included in a bulk p hase stability calculatio n, using 
Thermo-Calc software (Andersson et al. 2002). Eor a 
given set of Nc + 2 state variables [JPc being the num¬ 
ber of components) Thermo-Calc finds the global min¬ 
imum Gibbs free energy which lies on a plane tangent 
to the available phases’ Gibbs energy surfaces. A phase 
diagram representable as a 2d plot is then constructed 
from the set of tangent planes found by varying any two 
of the state variables. (Eor a pe dagogical refe rence to 
phase diagrams, see the book by Hillert (2007)). Since 
all the phases considered in this work are stoichiomet¬ 
ric crystal structures, there is a simplification in that 
the phases’ Gibbs energy surfaces are themselves points. 
Consequently, all phase regions in a phase diagram ob¬ 
tained by the above procedure must be 3-phase regions. 
If a structure appears in the equilibrium phase diagram 
at either P = 10^^ or 10^^ GPa, it is re-optimized at 
intermediate pressure decades to obtain the pressure de¬ 
pendence of the phase diagram. It is possible, though 
unlikely, that there are phases stable only over a nar¬ 
row band of pressure which are missed by this procedure 
(the full search scheme described in the last paragraph 


^ It is convenient to reinterpret XtalOpt and GULP’S internally- 
consistent (eV, A, GPa) unit system as (10^ eV, 10“^ A, 10^^ 
GPa), so that issues with numerical limits can be avoided. These 
codes were, of course, originally intended for Earth-condition ma¬ 
terials! Useful choices of the integer d include 2, 3 and 4. In this 
scheme, the relativity parameter appearing in Equations be¬ 
comes X = 1.1946484 X 10^^“^ , and the prefactor in Equation 

^becomes meC^/87r^Ag = 1.1239083 x 






























4 




Figure 1. Hartke plots for search cells Fe 03 C 2 (top) and 
Fe 404 C 4 (bottom). y{x) is the enthalpy of the lowest enthalpy 
structure found within the range of structure numbers: zero to x. 
The initial 20 seed structures (randomly generated by XtalOpt) are 
not included, thus the Hartke plot contains information only about 
structures which have been genetically-operated-on. For each plot, 
100 runs were made with identical parameters. In the case where 
all runs eventually found the same lowest enthalpy structure (top), 
the worst-best is the single run that took the longest to find it. 
In the case where not all runs found the same lowest enthalpy 
structure (bottom), the worst-best is the single run whose winning 
structure had the highest enthalpy. Best-best was the quickest to 
find the overall lowest enthalpy structure, and average-best is the 
average over all 100 runs. The Hartke lifetimes associated with the 
exponential fits to average-best (black dashed lines) are 35 (top) 
and 123 (bottom). The winning structure found in the Fe 03 C 2 
search appears in the bulk phase diagram, as the 77 phase. The 
Fe 404 C 4 search has the most degrees of freedom out of any search 
cell in our program; its winning structure does not appear in the 
phase diagram. 


was performed for the C-O-Fe system at several inter¬ 
mediate pressures and found no such “missed” phases, 
lending support to this approach). 

Five ternary systems of nuclei were selected for study: 
He-C-0, C-O-Ne, C-O-Fe, 0-Fe-Se, and Fe-As-Se. (The 
single exception to the search program described above 
has to do with He-C-0: our real space method converges 
much more slowly in this system due to the larger 
so the full 98 searches are carried out only at P = 10^^ 
GPa). The first two ternary systems are relevant to 
WDs, likely includin g those which are type l a super novae 
(SNIa) progenitors (Shen & Bildsten 2014). The third 
may also be relev ant to WDs having undergone a failed- 
detonation SNIa ( jJordan et al.||2Q12 ). The rationale for 
choosing the remaining two is that these partic ular nuclei 


are repres entative and/or prevalent among the Gupta et 


al. (2007) abundances near neutron drip. Incorporat¬ 
ing a lull list of abundances (~ 17 species) would be in¬ 
tractable for the type of calculation we have described. 
Moreover, we can take a lesson from earth-condition crys¬ 


tals, which typically have only 1, 2, or 3 elements (some¬ 
times 4). This appears to be due to general properties of 
crystal stability related to phase separation of complex 
unit cells: Basically, once a structure reaches a sufficient 
level of complexity that it can accommodate the spe¬ 
cial geometrical characteristics of its constituent atoms, 
it is disadvantageous to make the unit cell any larger 
(in the sense of adding more atoms), since that just re¬ 
duces the amount of favorable repetition possible with a 
given number of atoms. Exploring more ternary combi¬ 
nations is also likely to give diminishing returns in the 
prediction of new structures. Roughly speaking, with 
a smooth and spherically-symmetric interaction such as 
the Yukawa potential, there are only four qualitatively 
different ternary combinations: one big Z and two small 
Zs, two big and one small, all three mismatched, and all 
three similar. As long as the screening length regimes 
(characterized by /vs) are not too different, one ex¬ 
pects to see a continuity of structures formed by systems 
having similar relative Zs (or perhaps squared Zs), since 
there are no atom shell effects that come into play. An¬ 
other reason for studying the specific ternary systems 
mentioned above is that they cover at least three (ar¬ 
guably all four) of the qualitatively different combina¬ 
tions. 

3. BULK PHASE DIAGRAM RESULTS 

While pressure-invariance of the T = 0 phase diagram 
was anticipated in the extreme relativistic limit, it comes 
as some surprise that the pressure-independence persists 
well below this limit, nearly to the threshold for full pres¬ 
sure ionization. Figure shows that for all five ternary 
systems studied, no P-induced phase transitions were 
found above 10^^ GPa. For bee Fe, this pressure cor¬ 
responds to density 6.18 x 10^ g/cc and screening length 
kQ^/vs = 1.36, or about 75 percent of the saturation 
value. In general, screening length on the order of the lat¬ 
tice spacing appears to be a requisite for P-driven phase 
transitions. Further supporting this conclusion is the fact 
that the He-G-0 and G-O-Ne systems don’t undergo any 
pressure-induced transitions in the range lO^^-IO^^ GPa; 
within that range, screening lengths in these small Z sys¬ 
tems are significantly larger than one lattice spacing. 

Another finding is that combinations of nuclei with 
significantly mismatched Zs are much more conducive 
to efficient multicomponent packings than are systems 
where the Zs are fairly similar. For example, the Fe-As- 
Se system has an extremely simple low-pressure phase di¬ 
agram: at any composition, the microstructure consists 
simply of phase-separated bee crystallites. Multicompo¬ 
nent phases appear at high pressure, but they have the 
simple cesium chloride structure. In contrast, the G-O-Fe 
phase diagram is quite rich. Gombining one large Z and 
two smaller Zs results in a variety of binary and ternary 
crystal structures (enumerated in Table Q, all of which 
are more efficient (have a higher packing fraction) than 
phase-separated bee lattices. 

Both He-G-0 and 0-Fe-Se systems (two large, one 
small) feature all the same phases as G-O-Fe, except for 
the two ternary compounds which don’t appear. While 
a continuity of structures appearing between these sys¬ 
tems was anticipated, it is striking that at high pressures 
the two phase diagrams are identical. Glose similarity is 
also noted between the Fe-As-Se and G-O-Ne systems 
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Figure 2. T = 0 bulk phase diagrams for relativistic screened-Coulomb systems. Each pair of vertically-aligned diagrams corresponds to a 
specific ternary system of nuclei, while the two rows give the pressure dependence. Across the five ternary systems investigated, no pressure 
dependence was found in the range GPa, despite the fact that the screening length varies considerably over this range (from 

1.36rs to l.Slrs for bcc Fe). Following from our assumptions described in the main text, the microstructure in a given triangular region is 
a polycrystalline mixture of stoichiometric compounds (phases). All distinct phases are labeled with Greek letters and explained in Table 

m 


which both consist of three similar Zs. The outlier in 
this comparison is the nontrivial C-Ne binary structure, 
described in Tabled] These observations are consistent 
with the idea that it is the combination of relative Zs, 
and not of absolute Zs, that is important in determining 
the high-pressure phase diagram. 

In the pressure and screening length regimes appro¬ 
priate to this work (while P ranges from 10^^ to 10^^ 


GPa, /vs ranges from 1.13 to 1.81 for bcc Fe), there 
is a competition between close packing and next near¬ 
est neighbor interactions, which the closest packed struc¬ 
tures tend not to win. This is exemplified by bee’s favor- 
ability over fee, and the fact that only one of the equilib¬ 
rium phases found (magnesium diboride structure) also 
appears i n the phase diagram of de nsest binary sphere 
packings (Hopkins et al.|| 2011 | | 2012 ). The simplest mul- 
ticomponent crystals have structures that are also as¬ 
sumed by some ionic compounds under low pressure con¬ 
ditions, which may reflect the fact that ionic solids have 
a simple close-shell electronic structure (ionic solids also 
have strong +/— Coulomb interactions that are missing 
here). When a pair of Zs are not too dissimilar they 
usually form cesium chloride structure, e.g. OC, NeO, 
SeFe, SeAs and AsFe. When they are more dissimilar 
they tend to form magnesium diboride structure, e.g. 
OHe 2 , FeC 2 and Se 02 . Magnesium diboride is our first 
encounter with sub-cubic symmetry, which could give rise 


to transport anisotropy, elastic anisotropy, and other ef¬ 
fects such as a magnetic field coupling to the structure 
orientation. A quite prevalent but more complicated or¬ 
thorhombic structure occurs at chemical compositions 
04 He 4 , C 4 He 4 , Fe 404 , Fe 4 C 4 and 80404 , where these 
different instances can be interconverted by small adjust¬ 
ments of bond lengths and angles. A tetragonal struc¬ 
ture occurs at compositions C 4 He 2 and Fe 402 ; this is 
the second-highest density structure in the C-O-Fe sys¬ 
tem and could potentially drive oxygen to greater depths 
than it would otherwise go. The C-O-Fe system also fea¬ 
tures two ternary structures FeOC 4 and Fe 03 C 2 , both 
with hexagonal symmetry. FeOC 4 can be viewed as mag¬ 
nesium diboride structure, with the triangular magne¬ 
sium planes alternating between Fe and O compositions. 
Fe 03 C 2 consists of alternating layers of kagome O and 
honeycomb C, with Fe at the holes in the honeycomb 
layers. 

A general feature of the the ternary bulk phase dia¬ 
grams is that coexisting phases have mass density dif¬ 
ferences, due to a combination of neutron fraction and 
geometrical packing effects. These differences can be as 
large as ^ 10 percent of the total density and will result 
in stratification of phase domains in the presence of a 
gravitational field - the problem to which we now turn. 
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Table 1 

Selected compounds ajmearing in the C-O-Fe and C-O-Ne bulk phase diagrams, as indicated in Figure All numerical values given here 
correspond to P = 10^° GPa. Relative proton density is a measure of geometrical packing efficiency; r^tive baryon density includes the 
nongeometrical effect of neutron fractions. The reference phase for these relative densities is a-Fe, except in the case of 0-Ne2C4, for 
which the reference phase is a-Ne. Renderings have grey C, red O, green Fe, and violet Ne with sphere volume proportional to the nuclear 
charge Z. In the b and e renderings, bonds indicate Fe-Fe first nearest neighbors. If the space group is listed instead of a specific crystal 
structure, the unit formula gives the composition of the search cell in which the structure was found, not necessarily that of the primitive 
cell, pdb files of the structures (for all compositional instances) are included as supplementary materials in the online version. 


phase 

label 

unit 

formula 

crystallographic 
space group 

or structure 

relative 

proton 

density 

relative 

baryon 

densit3Q 

density relative 
to bulk, phase- 
separated bc(p| 

views along (or slightly oblique to) 
some high-symmetry directions 

a 

Fe 

bcc 

I 

I 

I 


a 

O 

bcc 

0.982 

0.912 

I 


a 

G 

bcc 

0.980 

0.910 

I 


d 

OG 

cesium-chloride 

0.981 

0.9II 

1.000001 


7 

FeG2 

magnesium-diboride 

0.994 

0.971 

1.000061 


5 

Fe 404 

Gmcm (orthorhombic) 

0.996 

0.979 

1.000040 

similar to (5-Fe4G4, see below 

»! » ! 

d 

Fe4G4 

Gmcm (orthorhombic) 

0.996 

0.983 

1.000056 






c 


FeOC 4 P6/mmm (hexagonal) 0.989 


0.950 1.000044 


• • • 

• • • 

••••••• 

•• •• 

• • • 





Fe 03 C 2 P6/mmm (hexagonal) 0.989 0.948 


1.000039 


••• • 
• _ • • 

• A 

• •• A • 

• • ••• 


• • • • 

•a«*a 
• • • • 

•a**a 
• • • • 

•a**a 


B Ne 2 C 4 Fd-3m (cubic) 


0.997 0.997 1.000021 


A A 9 9 9 9 

• • • • • • • 

• •(j|r ^ 9 9 99 

9 • • • • • • 

• • • • 

- JA- _ • • 

•a* •a* • • • 

A A 9999 


^Using 12c, 1^0, 20Ne, and ^^Fe 

^Defined as the sum of cell volumes after phase-separation into bulk bcc phases, divided by the original cell volume 


4. EQUILIBRIUM LAYERING GALGULATION 

Here we give an application of our high pressure crystal 
chemistry results to white dwarfs at a given fixed over¬ 
all composition. The equilibrium phase-layering diagram 
of a zero temperature WD is computed self-consistently, 
allowing for arbitrary numbers of components Nc and 
phases Np that can be formed from these components. 
The problem is decomposed into two parts: one part is 
a microscopic phase stability calculation which produces 
a function p{h) where p is density and h is enthalpy per 
unit mass, the other is a simple stellar structure calcu¬ 


lation which determines the pressure-radius dependence 
P(r). We iterate between these two parts. The former is 
inspired by a technique used among chemical engineers 
to study spe cies segregation in oil reservoirs, cf. [Esposito] 


et al. (20001) 

in the following, we will make use of the virial theorem 
for the gravitational potential energy W of a WD, given 
by 


nR 

W = —3 PAirr^dr. 

Jo 


( 7 ) 
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We begin by discretizing the star into Nl onion layers 
of uniform thickness Ar = R/N^. If Ar is small com¬ 
pared to the scale height of pressure Hp = —dr/dlogP, 
the layer may be treated as a bulk equilibrium sys¬ 
tem at constant pressure Pi, and one may work with the 
modified Helmholtz free energy 


Nl Np 

i=l a=l 


(8) 


For each term in the sum over layers, —SPiVi comes from 
the discrete version of Equation and another —PiVi 
cancels the corresponding quantity in the Gibbs free en¬ 
ergy of the layer, Uo^i Pi)- Here rid is the (un¬ 

known) molar amount of phase a present in layer i, and 
fia is the bulk chemical potential of phase a. (The phase 
index a is not to be confused with the bcc structure, as 
in Table [^. We have thus avoided the complication of 
introducing a gravitational term into the chemical po¬ 
tentials, including it instead at the level of the layers. 
This comes at the cost of supplying a pressure function 
P(r) consistent with hydrostatic equilibrium, implicit in 
Equations [^&[^ Let’s assume that we have such a pres¬ 
sure function. (Eor an initial guess, we will take P(r) 
from a n = 3 poly trope.) Now fix a set of layer pressures 
Pi = P{iAr). The problem of minimizing P* has been 
reduced to the problem of minimizing the linear objective 
function rid /^a{T, Pi), subject to 2Nl P Nc - I 

constraints 


~V,pPa{T,Pi) 

0 = y][(l -a;A)sAa 

i,a 

etc. for xp - - - xnc^ 


for i = 1... Nl, 
^A^Ba ^A^Ca] '^ai, 


( 9 ) 


( 10 ) 


0 > 


i:[ 


Vi+i 


'^ai 


rrio 


fori = l...AL-l,(ll) 


each of which is also linear in the rid- The first set of 
constraints ensures the volume filling fraction is equal to 
1 for each layer, where the molar mass rria and density 
pa(P^Pj) assumed to be known for each phase. Equa¬ 
tions [M constrain the global mole fractions xa- - - 
whereTTor example, SAa specifies the number of H-type 
nuclei per formula unit of the a phase. The reason 
for constraining mole fractions rather than component 
masses is that the latter tends to cause infeasibility prob¬ 
lems for the Simplex solver. Einally, there is a set of in¬ 
equality constraints which guarantee reality of the inter¬ 
layer Brunt-Vaisala frequencies uj = ^/—{g/p){dp/dr). 
Thus, there is built-in stability against convective over¬ 
turn of adjacent layers, but note that it is still possible to 
have unstably stratified material within a layer. As noted 
above, this problem is straightforwardly solved using the 
Simplex method. Eor number of variables NlNp ^ 10^- 
10^, we use the high-performace lp_solve routines. 

So far we have considered the case of an isothermal 
WD. Eor the special case T = 0, the Simplex solu¬ 
tion provides the layer enthalpy per unit mass hi = 

^ud layer density pi = 
In Certain cases, one can interpolate 


to obtain a smooth function p{h), which can be com¬ 
bined wit h the enthalpy- transformed stellar structure 
equations (Lindblom|l992 ) ^ In the nonrelativistic limit, 
these read 


dP 

dm —Airr^p 

(12) 

(13) 

dh Gm ’ 

dr —r^ 

(14) 

dh Gm 


The reason for using the enthalpy transformation is 
twofold. Eirst, the total mass M, which we were not 
able to constrain in the Simplex calculation, now enters 
as a boundary condition. Second, if we simply used the 
layer masses Mi = Uai rua along with the discretized 
equations of mass continuity and hydrostatic equilibrium 
(the usual stellar structure equations) to update the Pi, 
no information about the microscopic energy scale would 
carry over. In other words, we could multiply all the 
chemical potentials by 2 and get the same Pi. Equations 
p!2}|T^ are to be integrated inward from the boundary 
conditions P{0) = 0, m(0) = M, and r(0) = R. Un¬ 
fortunately we have no apriori knowledge of the radius 
R that is consistent with M and p{h) in the sense of 
the Oppenheimer-Volkoff map (to use the language of 
Lindblom). We therefore have to complete the mapping 
p{h) \-A {M,R). A simple way to accomplish this is by 
“aiming” the boundary condition r(0) until the integra¬ 
tion yields the physically-correct behavior dP/dr = 0 as 
r ^ 0. Approaching R from below, the solutions are 
smooth and well-behaved except at r = 0 due to a sin¬ 
gularity in Equation M A change of variable u = r‘^ 
removes this singularity, but we find no particular ad¬ 
vantage to working with the resulting transformed equa¬ 
tions. Approaching R from above generates sign changes 
and the solutions are generally chaotic. The qualitatively 
different behaviors in these two regimes can be exploited 
to obtain R to arbitrarily high precision. In practice, we 
minimize dP/dr at a fixed, small fraction of the starting 
boundary condition r(0) (but see the next paragraph for 
discussion of a special case). In this process of “com¬ 
pleting the map,” an updated pressure function P(r) is 
obtained at no extra cost. Layer pressures are reassigned 
and input to the Simplex calculation, and the process 
iterated. One choice of convergence criterion is that suc¬ 
cessive iterations produce stellar radii which are the same 
to within a tolerance of 10 “^Pq - typically this criterion 
is met within just a few iterations. Another choice is 
that radial positions and thicknesses of phase strata (as 
fractions of R) are static to within the resolution set by 


^ An issue can arise near a density discontinuity, where hi and pi 
obtained by the above procedure describe a function h{p) which is 
non-monotonic. The stellar structure calculation cannot then m ake 
use of the enthalpy transformation, because the sign of Equation |14| 
is incorrect in the vicinity of the interface. A density discontinuity 
occurs as a consequence of mismatched Zs (compare bcc C and 
O in Table but is made much more severe when there is also 
a mismatch m neutron fraction (compare bcc O and Fe in Table 
[n. Fortunately, for typical white dwarf compositions, the neutron 
maction is continuous (or nearly so) across phase boundaries and 
the issue of non-monotonicity is avoided by choosing a suitably 
large layer thickness - on the order of i^/200 for He-C-O and C-O- 
Ne compositions. 













the number of simulation layers - typically this occurs 
after just one iteration. 

The main type of numerical error incurred is of the fol¬ 
lowing nature. In the first iteration, integration of Equa¬ 
tions 12 -14 never proceeds past the point for which we 
have tabulated pi data available to interpolate within. 
This is just a consequence of having used the polytrope 
initial guess. In subsequent iterations, however, we are 
sometimes forced to make a choice: carry out the inte¬ 
gration past the highest tabulated replacing interpo¬ 
lations with extrapolations, or simply terminate the inte¬ 
gration when interpolations become impossible, using the 
current value of dP/dr in the aiming procedure discussed 
above. We choose the second option. The miminization 
problem within the aiming procedure is, in these cases, 
somewhat ill-defined, tending to generate some numeri¬ 
cal noise which is expressed in the layering diagram near 
r/R = 0. For this reason we present layering diagrams as 
they appear after the first iteration, noting that changes 
to the layering diagram are already imperceptible by the 
second iteration, save for an increase in the level of this 
numerical noise. 


5. EQUILIBRIUM LAYERING RESULTS 

The previous section described a method of determin¬ 
ing the radial positions and amounts of phase strata 
in a white dwarf with fixed mass and overall composi¬ 
tion, in particular, strata composed of the new multi- 
component crystal structures. Using this method, we 
computed the T = 0 equilibrium phase layering dia¬ 
grams and radius-composition dependence of 1.0 solar 
mass, ^He-^^C-^^0 and white dwarfs. For 

each composition, an initial guess corresponding to the 
polytrope P = 3.8 x c.g.s. and stellar radius 

P = 7.5 X 1 O“^P 0 was used, although the method ap¬ 
pears to converge to the same result if these starting val¬ 
ues are adjusted within reasonable limits. Figureshows 
the result of the calculation. Evidently pure bee phases 
make up the majority of the stellar interior, despite the 
multicomponent structures being more efficiently packed. 
Since multicomponent phases tend to show up at inter¬ 
faces, we refer to them as “interphases.” In the He-C-0 
star, for example, J-C 4 He 4 and e-C 4 He 2 interphases are 
formed between a-C and a-He, while /3-OC appears at 
the low density part of the C-0 boundary. For composi¬ 
tions near xpe = 1 ? fhe thinness of the carbon shell allows 
0-He interphases to form, namely 7 -OHe 2 and ( 5 - 04 He 4 . 
Compared to sharp bcc-bcc interfaces, interphases offer 
free energy savings due to optimized crystal packing den¬ 
sity, nearest and next-nearest neighbor interactions, etc. 
arising from the extra compositional degrees of freedom. 
Interphase thinness relative to bee strata can be under¬ 
stood from the gravitational contribution to the free en¬ 
ergy having a tendency to “pull apart” the different Z 
components of the multicomponent phases. Competition 
between these two energy scales apparently causes inter¬ 
phases to become slightly thicker with depth. Consider, 
for example, ( 5 -C 4 He 4 in the diagram with xqIxq = 1 . 
At XHe = 0.5, only one simulation layer (out of 200 ) is 
completely filled with this compound, while an adjacent 
layer contains a mixture of ^-C 4 He 4 and a-He. This in¬ 
terphase gradually thickens with xpe by xpe = 0.95, 
13 simulation layers are completely filled, another two 
are partially filled, and (i-C 4 lIe 4 has squeezed out a-C 


from the layering diagram. 

An unexpected but apparently generic feature of the 
radius-composition curves in Figure is the existence of 
a shallow minimum of the WD radius at an impure com¬ 
position. Even as the level of numerical noise increases 
with further iterations, this minimum clearly persists. 
The cusps near xpe = 0.8 and x^q = 0.3 may be a nu¬ 
merical effect rather than a physical one, however, as 
they tend to smooth out upon further iterations. 

6 . NONEQUILIBRIUM LAYERING RESULTS 

A simple modification of the equilibrium layering cal¬ 
culation enables a quasi-static settling calculation. If 
the settling species is X, an additional set of linear con¬ 
straints 

0 — ^ ^ ^Xa'^ai i 
a 

enforces the minimum depth imin at which X can appear. 
This minimum allowed depth can then be incrementally 
stepped down. We carry out a test of the method by set¬ 
tling O.OOMq of O on a O.OIMq He-C white dwarf, and 
0.1 Mq of Ne on a 0.9 Mq C-0 white dwarf. Overall com¬ 
positions are fixed at xpe = 0.95 with xq = xq = 0.025, 
and = 0.07 with xq = 2xq = 0.62, respectively; the 
final settled-out states are given by Figure ^ (Admit¬ 
tedly, these are not particularly realistic settnng scenar¬ 
ios but they serve as interesting test cases, forcing the 
presence of an interface between the highest and lowest 
Zs which otherwise doesn’t happen in equilibrium). Note 
that if the starting value for imin is too near the surface, 
there is no feasible solution that can accommodate the 
settling mass. For this reason we restrict our study to the 
second half of the settling problem: imin = ..., 0 . 

Results of the quasi-static settling calculation are plot¬ 
ted in Figure In both settling scenarios, the out-of- 
equilibrium star contains one or more phases that do not 
appear in the final, equilibrium stacking sequence. One 
function of these extra phases is to serve as transient host 
structures for the settling species: X 04 He 4 and 7 -OHe 2 
are hosts for settling oxygen, and 6 >-Ne 2 C 4 is a host for 
settling neon. The phase settling diagram for the He-0- 
C star is fairly non-trivial, with as many as seven distinct 
strata near imin/^L ~ 0.25. A minimum in the stellar 
radius appears around this point (as it also does in the O- 
C-Ne settling calculation) indicating that the lowest en¬ 
thalpy star is not the most compact star. This minimum 
appears to track the size of the e-C 4 He 2 core, possibly 
also the thickness of the X 04 He 4 interphase. However, 
the minimum persists when the calculation is repeated 
excluding first one and then the other of these phases. 
These result hints at the prospects for new phenomena 
that are enabled by compositional and structural hetero¬ 
geneity that takes advantage of the additional “chemical” 
degrees of freedom afforded by multinary phases. 

In the first settling calculation, both He and C 
must eventually find their way through the sinking O- 
containing layer. It is interesting that, with the exception 
of a single point near iminj^L =0.3 where all the oxygen 
is bound up in X 04 He 4 , there is no continuous migra¬ 
tion pathway (in the sense of stoichiometric compounds) 
assuming any deviations from spherical symmetry are 
sufficiently weak to maintain contiguity of the layering 
sequence. Several binary phases are available to pro- 
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0.75 


r/R 0.5 


0.25 


0 0.25 0.5 0.75 1.0 

XHe (xo/xc=1) 


0 0.25 0.5 0.75 1.0 

XHe (xo/xc=2) 





0 0.25 0.5 0.75 1.0 

XHe (xo/xc=4) 



XNe (Xo/Xc=0.5) 



XNe (Xo/Xc=1) XNe (Xo/Xc=2) 


Figure 3. Equilibrium phase layering diagrams for I.OMq, He-C-O white dwarfs (top row) and I.OMq, C-O-Ne white dwarfs (bottom 
row). Carbon-oxygen ratio is held at a fixed value in each panel, and mole fraction He or Ne makes up the balance. The discrete color map 
indicates the stable phases after one iteration of the procedure described in the main text. Further iterations produce no perceptible change 
in the layering diagram, save for an increase in the level of numerical noise. In cases where a given layer contains a two phase mixture of 
bcc and a more complicated structure (these are the only types of mixture to occur), the color map indicates the more complicated phase. 
White crosses give the stellar radius (right y-axis) as a function of composition, also after one iteration of the method. 200 layers were used 
in the computation of these diagrams. 










10 



imin/NL (time x settling rate) 



Figure 4. Quasi-static nonequilibrium phase diagram for settling 
0.09Mq of O on a 0.91Mq He-C white dwarf (top) and settling 
O.IMq of Ne on a 0.9Mq C-O white dwarf (bottom). The ar-axis 
can be regarded as time remaining until equilibrium, multiplied by 
the settling rate. White crosses give the evolution of the stellar 
radius during the settling process. The top (bottom) diagram was 
computed using 180(200) layers. Other details of the plots are the 
same as in Figure 


vide such a pathway, but the system does not use them 
to this advantage - for example, /3-OC is never formed. 
Consequently, in the final stages of settling, carbon has 
to diffuse through an oxygen barrier having thickness 
t ^ 10^ km. The associated timescale can be estimated 
as r ^ where Dq = is t he diffusion 


coefficien t for a one-compo nent plasma (see Ha nsen et 
al. ( 1975| ) and more recently, Hughto et al | ( ^iu[ )). Here 


ftp = ZnelMy^‘^ is the ion plasma frequency and 
r = Z‘^e‘^/{vsksT) is the Coulomb coupling parameter. 


which we take to be the melting point value: T^ = 175. 
Putting in the numbers gives r ^ 10^^ yrs. While this 
is an oversimplified analysis, it does suggests that com¬ 
pact object phase strata may be far from the equilib¬ 
rium stacking sequence. It is particularly intriguing to 
consider whether strong “chemical” deviations from the 
equilibrium phase stacking sequence could accumulate 
excess free energy that is eventually liberated in ener¬ 
getic (and thus observable) events. Our analysis here 
is a first step towards developing a framework towards 
consideration of such possibilities. 


7. FINITE TEMPERATURE 

Here we give a brief, qualitative discussion of some ef¬ 
fects that will become important at finite temperatures; 
a detailed analysis is a subject for future work. At finite 
T, the chemical potentials fia{T,Pi) must be modified 
to account for smearing of the Fermi surface, associated 
smearing of the electron capture layers which will give 
an adjustment in composition, phonons, and if a de¬ 
notes an alloy or solution phase, mixing entropy. Since 
the entropic part of these contributions does not enter 
into the enthalpy, it is not possible to self-consistently 
include thermal effects in our equilibrium phase layer¬ 
ing method, which relies on the enthalpy transformation. 
However, thermal effects could be included post hoc. For 
example, one could compute the phonon free energy of 
an interphase crystal such as (5-C4He4 along with that of 
the equivalent phase-separated a-C and a-He crystals, 
add this quantity to the T = 0 free energy, and pre¬ 
dict whether the interphase tends to thicken or thin at 
finite T. A rough estimate based on phases’ bulk moduli 
{K = —VdP/dV) suggests the thermal (phonon) correc¬ 
tion to the free energy will tend to increase the stability 
of the soft outer phase strata relative to the stiff innner 
strata, likely shifting the interphases slightly towards the 
stellar center. There is also a possibility for additional, 
particularly soft interphases to appear in the stacking 
sequence, if entropic terms are large enough to affect 
the phase competition that winnows the phases of Ta¬ 
ble [l] down to the phase layering in Figure A phonon 
calculation would also yield the elastic tensor and help 
characterize the degree of elastic anisotropy as a func¬ 
tion of depth, as well as give a prediction for the relative 
melting temperatures via the Lindemann parameter (ra¬ 
tio of RMS nucleus displacement to equilibrium lattice 
spacing). Finally, we note that the simple mechanical 
stability criterion used here shou ld be replaced with the 
Ledoux criterion at finite T. See Reisenegger (2001) for 
an application to multicomponent neutron stars. Again, 
it is difficult to see how this more sophisticated criterion 
can be built in to the current method, but it could at 
least be checked after including thermal effects in the 
manner described. 
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